RIF and PMB combination

Background

We looked for combination therapies of neglected and misused antibiotics. Rifampin and polymyxin B combination stood out as a promising approach against scpectrum of clinical isolates. This repository contains all the data and code to reproduce the analysis. For details and full captions of figures, consult the publication.

Combinations

We studied the potential of the rifampicin-polymyxin B combination against intra- and extracellular forms of bacteria: three P. aeruginosa strains, two clinical isolates of A. baumannii, E. coli, and K. pneumoniae.

dat = loadDdiData('input/dat/rds/ddi_all_strains.rds')

# Response surface analysis
# rs = getRS(dat, 'rs')
# rs = readRDS('input/dat/tmp/rs.rds')
# 
# open3d()
# lapply(names(rs), sav3D)

#getDdiCounts(readRDS('input/dat/tmp/rs.rds')) %>% 
#  group_by(Strain) %>% 
#  ggplot(aes(Cond, Syn)) +
#  geom_col(fill = 'grey70', col = 1) +
#  facet_wrap(~Name, nrow = 3, dir = 'v') +
#  labs(y = '# of synergistic concentrations', x = '') +
#  theme(strip.text.x = element_text(size = 12)) 
#ggsave(filename='output/fig/SynergyCountAcid.svg', bg='white', h=7.5, w=7.5)

Dose Response And pH

theme_set(theme_pdf())

sub = subset(dat, strain == 'ATCC27853' & (d1 == 0 | d2 == 0))

p1 = function() pltDR(effect ~ c1, subset(sub, c2==0), xlab = 'RIF, mg/L', zero=30)
p2 = function() pltDR(effect ~ c2, subset(sub, c1==0), col='orange', xlab = 'PMB, mg/L')

# pH
foo = fread("input/dat/raw/growth_and_pH.tsv")

foo = foo[time_h >= 0 & overnight_medium == 'pH7.4', .(
  od = gmu(od620),
  od.ci = gci(od620),
  ph = mean(pH),
  ph.ci = ci(pH)
  ), .(growth_medium, time_h)]

limit_of_detection = 0.03
foo[, od := replace(od, od < limit_of_detection, limit_of_detection)]

p0 = foo %>% 
  mutate(
    growth_medium = relevel(factor(growth_medium), 'pH7.4'),
    time_h = ifelse(time_h > 20, 20, time_h),
    NULL
    ) %>% 
  ggplot(aes(time_h, od, shape = growth_medium)) +
  scale_x_continuous(breaks = c(0, 5, 10, 20), labels = c(seq(0, 10, 5), 24), limits = c(0, 20)) +
  scale_shape_manual(name = '', values = c(19, 22, 15), labels = c('pH 7.4', 'pH 5.5',
      'pH 5.5 buffered')) +
  theme(legend.position = 'top')

p3 = p0 + 
  scale_y_continuous(trans='log2', breaks = c(0.01, 0.03, 0.1, 0.3, 1, 3,
      10)) + 
  geom_point(size=.point_size, col='grey20') +
  geom_pointrange(aes(ymin = od*od.ci, ymax = od/od.ci)) +
  geom_line(col='grey20') +
  coord_cartesian(ylim = c(0.01, 10)) +
  labs(x = 'Time, h', y = expression('OD'[620]))

p4 = p0 %+% aes(y = ph) + 
  scale_y_continuous(position = 'right') +
  geom_point(size=.point_size, col=cbPalette[3]) +
  geom_line(col=cbPalette[3]) +
  geom_pointrange(aes(ymin = ph-ph.ci, ymax = ph + ph.ci), col=cbPalette[3]) +
  theme(legend.position = 'none',
    axis.title.y.right=element_text(color = cbPalette[3]),
    axis.text.y.right=element_text(color = cbPalette[3]),
    axis.ticks.y.right=element_line(color = cbPalette[3])
  ) +
  labs(x = '', y = 'pH')

foo = fread('input/dat/raw/charcoal.csv')
labs = c('RIF, mg/L', 'PMB, mg/L', 
  rep(expression(log[10]*'('*over('CFU'['24h'], 'CFU'['0h '])*')'), 2)
)

tmpFun = function (i, df) {
  par(mar=c(3.1, 6.1, 0, 0.7), mgp = c(1.7, 0.5, 0), tcl=-0.2)
  if(i == 1) bar = subset(df, ab == 'rif')
  else bar = subset(df, ab == 'pmb' | (ab == 'rif' & conc == 0))

  # fit char vs no-char but share Emin and EC50 (argument pmodels)
  fit = drm(response~conc, data=bar, CURVE, 
    pmodels=data.frame(CURVE,1,1),  
    fct=LL2.4( fixed=c(1, NA, NA, NA))) 
  
  plot(fit, type="average", cex.axis=1.25, cex.lab=1.25, lwd=2, 
    col='grey30', pch = c(1, 19), lty = 1, asp=1,
    xlab = labs[i], ylab = labs[i+2], broken=T,
    ylim=c(-4,4), legendText=c("control", "charcoal")
  )
  plot(fit, type="confidence", add=T, legend=F, 
    col=c('grey70', 'grey10')
  )
}

p5 = function() tmpFun(1, foo)
p6 = function() tmpFun(2, foo)

aligned = align_plots(p3, p4, align = 'hv', axis = 'tblr')

plot_grid(
  ggdraw(p1), ggdraw(p2), ggdraw(aligned[[1]]) + draw_plot(aligned[[2]]),
  plot_grid(ggdraw(p5), ggdraw(p6), nrow=2),
  labels = 'AUTO', align = 'v'
)

# ggsave('output/fig/SFig_mth_dr_pH_charcoal.svg', width=10, height=9)

PD parameters

# schematic for PD pars
tmp = data.table(
  #Conc      = c(0, 0.1, 0.3, 1, 3, 10, 30, 100),
  Conc      = c(0, 0.01, 0.03, 0.1, 0.3, 1, 3, 10),
  Response  = c(3, 3, 2.8, 1.5, -0.5, -2.5, -2.8, -3.2)
)

tmp[, Response := Response + 1]

# svg('output/fig/SFig_PDpar.svg', width = 3.5, height = 3.5)
#oldpar = par()
par(tcl = -0.3, mar=c(4.1, 4.1, 1.5, 0.5), mgp=c(1.5, 0.75, 0))

plot(
  (mod = drm(Response ~ Conc, data = tmp, 
      fct=LL2.4(names=c("h", "m", "b", "e"), fixed=c(1, NA, NA, NA)),
      lowerl = c(-2.5, NA, NA),
      robust = 'median'
  )),
    broken = T, pch = 19, 
    col = 'grey40',
    ylim = c(-2.5, 4.7), xlim = c(0, 50),
    xlab = c('Conc, mg/L'),
    ylab = expression(log[10]*'('*over('CFU'['24h'], 'CFU'['0h '])*')'),
    asp = 1
)

x1 = ED(mod, 0, type = 'absolute', logBase = exp(1))[1]

## 
## Estimated effective doses
## 
##        Estimate Std. Error
## e:1:0 0.4249298  0.0092374

x2 = ED(mod, 50, logBase = exp(1))[1]

## 
## Estimated effective doses
## 
##        Estimate Std. Error
## e:1:50 0.233915   0.005085

y2 = predict(mod, data.frame(Conc = x2))

segments(x1, 0, x1, -100, lty = 2)
segments(0.0001, 0, x1, 0, lty = 2)
segments(x2, y2, x2, -100, lty = 2)
segments(0.0001, y2, x2, y2, lty = 2)

abline(h = coef(mod)[2], lty = 2)
abline(h = coef(mod)[1], lty = 2)

text(x1+0.3, y = 0 + 0.1, expression(C[s]))
text(x2+0.3, y = y2 + 0.2, expression(EC[50]))
text(x = 0.002, coef(mod)[2] + 0.3, expression(E[min]))
text(x = 0.002, coef(mod)[1] + 0.4, expression(E[max]))

# dev.off()

# add Tn screen's PA14 (Liberati), and pmbR strains PA947 and PA1292 
foo = rbind(
  # pa14
  fread('input/dat/raw/pa14_mut_pmb_rif_cfu.csv')[mut == 0] %>% 
    .[, `:=` (genus   = 'Pseudomonas',
              strain  = 'PA14 Liberati',
              mut     = NULL)],
  # pmbR
  fread('input/dat/raw/pmb_resistant_isolates.csv') %>% 
    .[, genus := 'Pseudomonas']
  ) %>% 
  .[, `:=` (
    c1 = d1*mic1,
    c2 = d2*mic2)]

# transform
dat_fit = rbind(foo, dat, fill = T)[c1 == 0 | c2 == 0,
  ][, dose := ifelse(d1 == 0, c2, c1)
  ][, ab   := ifelse(d1 == 0, 2, 1)
  ][, .(data = list(.SD)), .(strain, cond)]

# fit four parameter logistic regression
dat_fit[, fit := map(data, ~drm(effect ~ dose, ab, data = .x, 
  fct=LL2.4(names=c("h", "m", "b", "e"), fixed=c(1, NA, NA, NA)),
  pmodels = data.frame(ab, 1, ab),
  lowerl = c(-4.5, NA, NA)))]

# get parameters
dat_fit[, pars  := map(fit, ~getPDPar(.x))]

## 
## Estimated effective doses
## 
##        Estimate     Lower     Upper
## e:2:0    1.0936    0.7459    1.6033
## e:1:0 3017.5829 1432.1220 6358.2618
## 
## Estimated effective doses
## 
##        Estimate     Lower     Upper
## e:1:0 199.80084 143.78167 277.64579
## e:2:0   1.09877   0.85522   1.41167
## 
## Estimated effective doses
## 
##         Estimate      Lower      Upper
## e:1:0  4195.4112   844.9468 20831.4600
## e:2:0     1.4129     1.1031     1.8096
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  136.099  94.119 196.802
## e:2:0  108.833  77.678 152.485
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  116.627  78.836 172.534
## e:2:0  100.888  71.106 143.146
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  2.19824 1.73299 2.78839
## e:2:0  1.05512 0.74755 1.48923
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:2:0  0.54794 0.41895 0.71664
## e:1:0  2.54461 2.01152 3.21898
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:2:0  1.23558  0.90005  1.69619
## e:1:0 12.27513  9.50191 15.85774
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:2:0  0.72807  0.50460  1.05050
## e:1:0  8.19634  6.45160 10.41292
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:2:0   8.8108  5.2898 14.6755
## e:1:0  29.9697 20.8806 43.0152
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:2:0   2.9223  1.9854  4.3015
## e:1:0  25.5587 19.0791 34.2388
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:2:0  1.08515  0.83477  1.41063
## e:1:0 17.68410 13.82126 22.62653
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:2:0   3.0869  1.4710  6.4777
## e:1:0   9.1434  7.0217 11.9061
## 
## Estimated effective doses
## 
##       Estimate  Lower  Upper
## e:2:0   2.7579 1.6050 4.7391
## e:1:0   6.2385 4.5268 8.5975
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  0.71158 0.43400 1.16668
## e:2:0  4.78692 3.03655 7.54625
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0  0.61538  0.37753  1.00308
## e:2:0 11.63563  7.65409 17.68830
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 10.31467  8.03876 13.23493
## e:2:0  0.52524  0.40478  0.68156
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 21.90092 17.02990 28.16518
## e:2:0  0.58875  0.45974  0.75395
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0   34.814  26.015  46.590
## e:2:0  212.935 151.138 300.000
## 
## Estimated effective doses
## 
##       Estimate  Lower  Upper
## e:1:0   55.495 42.649 72.210
## e:2:0   24.544 18.867 31.930
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  22.7675 16.2560 31.8872
## e:2:0   2.0841  1.4483  2.9992
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 18.17315 12.63799 26.13261
## e:2:0  0.36616  0.23469  0.57127
## 
## Estimated effective doses
## 
##        Estimate     Lower     Upper
## e:1:0  6.273396  3.311455 11.884652
## e:2:0  0.071390  0.039096  0.130359
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  3.89931 2.70443 5.62209
## e:2:0  0.50220 0.38022 0.66331
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  5.68227 3.61600 8.92927
## e:2:0  1.27999 0.91009 1.80022
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  15.6419 11.5468 21.1894
## e:2:0   1.1580  0.8129  1.6495
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 34.12972 26.34817 44.20944
## e:2:0  0.56378  0.45798  0.69403
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  63.2718 45.7219 87.5581
## e:2:0   6.4509  4.8977  8.4967
## 
## Estimated effective doses
## 
##        Estimate     Lower     Upper
## e:1:0  78.80048  58.25397 106.59387
## e:2:0   1.15853   0.89306   1.50291
## 
## Estimated effective doses
## 
##       Estimate   Lower   Upper
## e:1:0  19.0752 14.1575 25.7011
## e:2:0   3.9380  2.8444  5.4521
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:1:0 38.15049 27.24067 53.42967
## e:2:0  0.92200  0.66327  1.28166
## 
## Estimated effective doses
## 
##       Estimate    Lower    Upper
## e:2:0  0.60558  0.38290  0.95774
## e:1:0 27.16329 14.77923 49.92439

dat_fit = dat_fit[
  dat_fit[, data[[1]][1, .(genus, mic1, mic2)], .(strain, cond)], 
  on = c('strain', 'cond')]

# reformat
dat_fit = unnest(dat_fit, pars) %>% as.data.table
dat_fit = getPdReadyForPlotting(dat_fit, dat)

p1 = dat_fit[!grepl('C', par)] %>% 
  .[lwr < -10, c('lwr', 'upr', 'value') := NA] %>% 
  ggplot(aes(value, strain, color = genus)) +
    geom_vline(xintercept = 0, col = 'white') +
    geom_pointrange(aes(xmin = lwr, xmax = upr), size = 1) +
    geom_text(aes(label = round(signif(value, 2), 1)), col = 1, size = 2.5) +
    facet_grid(cond ~ par, labeller = label_parsed, 
      scales = 'free', space = 'free_y') +
    geom_blank(aes(x = min)) +
    geom_blank(aes(x = max)) +
    scale_color_manual(values = cbPalette[4:1], guide = guide_legend(
        label.theme = element_text(angle = 0, face = "italic"), reverse = T)) +
    scale_fill_manual(values = cbPalette[4:1]) +
    labs(y = '', color = '') +
    theme_html() +
    theme(legend.position = 'top')

p2 = p1 %+%
  mutate(dat_fit[grepl('EC', par)], 
    value = ifelse(value > 512, 512, value), 
    mic = ifelse(mic > 512, 512, mic)) +
  geom_point(aes(mic, col = genus), size = 4.5, pch = 0) +
  geom_text(aes(x = mic, label = round(mic, 1)), size = 2.5) +
  scale_x_continuous(trans = 'log10', limits = c(0.03, 512), 
    labels = function(x) ifelse(x == 0, "0", x)) +
  facet_grid(cond ~ par, labeller = label_parsed, 
    scales = 'free', space = 'free_y') +
  xlab('Conc, mg/L') +
  theme(legend.position = 'none')

plot_grid(p1 + xlab(expression(log[10]*'('*over('CFU'['24h'], 'CFU'['0h '])*')')), 
  p2, rel_widths = c(3, 2), align = 'hv', axis = 'tb')

# ggsave('output/fig/SFig_mth_pd.svg', width=10, height=6)

Time-kill

dtk = fread('input/dat/raw/time-kill.csv') %>% 
  remNonCharcoal() %>% 
  addNewVariablesToTimeKill()

# pltTimeKill(dtk)
#ggsave('output/fig/SFig_TimeKill.svg', width=15, height=5)

# # Print out the concentrations used
# unique(dtk[time_h > 0], by = c('broth', 'group')) %>% 
#   .[, .(broth, group, d_1, d_2)] %>% 
#   .[order(group)]

Reverse Genetics Screen

To account for the synergy in molecular terms—beyond a nonspecific increase in membrane permeability by polymyxin B—we turned to chemical genetics (Brochado and Typas, 2013). Working with ordered PA14 transposon library (Liberati et al., 2006), we derived a growth measure for monotherapies and combinations using colony opacity (Kritikos et al., 2017). To account for plate-to-plate variation, the opacity was multiplicatively corrected. This results in zero-centering of the Bliss scores, which were derived next. The significance of difference from zero Bliss score, for any mutant, was estimated by a T-test (5 biological replicates) and corrected for multiple testing (Benjamini-Hochberg).

tmp = chem.gen.res %>% 
  mutate(t.test.q.value = ifelse(t.test.q.value < 0.001, 0.001,
      t.test.q.value)) %>% 
  mutate(showname = ifelse(
      (t.test.q.value < 0.05 & (gene.name != '' | gene.name.to.show %in% c('66480', '26590', '43270', '02150'))) | t.test.q.value < 0.01 ,
      T, F)) 

# tmp[t.test.q.value < 0.05 | t.test.q.value.gene < 0.05, .N, media]
# #      media  N
# # 1:      LB 65
# # 2: LBpH5.5 39

# # only 3 are shared between the media
# foo = tmp[media == 'LB' & (t.test.q.value < 0.05 | t.test.q.value.gene < 0.05), Tn.mutant.id] 
# bar = tmp[media == 'LBpH5.5' & (t.test.q.value < 0.05 | t.test.q.value.gene < 0.05), Tn.mutant.id] 
# tmp[Tn.mutant.id %in% intersect(foo, bar)]


getLocus = function (dat, condition) {
  unique(dat[eval(parse(text=condition))], by = c('locus'))[, locus] 
}

by_gene       = getLocus(tmp, 't.test.q.value.gene < 0.05')
by_mut        = getLocus(tmp, 't.test.q.value < 0.05')
by_both       = getLocus(tmp, 't.test.q.value < 0.05 | t.test.q.value.gene < 0.05')
by_gene_only  = setdiff(by_gene, by_mut)

# # 70 hits
# length(by_both)
# # 53 from mutants
# intersect(by_mut, by_both) %>% length()
# # 17 additional from genes
# setdiff(by_gene, by_mut) %>% length() 
# # 31 shared between mutants and genes 
# intersect(by_gene, by_mut) %>% length()

Of the 70 hits, 53 are from mutant level of analysis and 17 come from gene level of analysis. 31 are shared between mutants and genes i.e. hits regardless if mutant or gene level of analysis is used.

(p0 = filter(tmp, media == 'LB') %>%  
  {
  ggplot(., aes(epsilon.median.mutant.condition, -log10(t.test.q.value), col = t.test.q.value < 0.05)) +
  geom_point(size = 2.5) +
  scale_color_manual(values = c('grey70', 'grey20')) +
  ggrepel::geom_text_repel(data= . %>% .[(showname)], 
    aes(label=gene.name.to.show), size = 6,
    max.overlaps = Inf, box.padding=0.75, seed=3) +
  lims(x = c(-0.8, 0.6), y = c(0, 3)) +
  labs(x = 'Shift in Bliss score', y = expression('log'[10]*'(p-value)')) +
  theme(aspect.ratio = 1, legend.position = 'none')

})

# ggsave('output/fig/ChemGen_LB.svg', width=5, height=5)

p0 %+% filter(tmp, media == 'LBpH5.5')

# ggsave('output/fig/ChemGen_LBpH5.5.svg', width=5, height=5)

LB	LB pH 5.5

getGenes = function (dat) {
  out = list()
  for(i in unique(dat$media))
    out[[i]] = unique(dat[t.test.q.value < 0.05 & media == i], by='locus') %>%
      .[, (gene.name.to.show)]
  return(out)
}

foo = getGenes(chem.gen.res)
#vtable = gplots::venn(foo)

GO terms

With the following, we bring some biological knowledge into the analysis. This will get us at the level of processes/compartments as opposed to individual genes. We will focus on LB beacuse we have more data from there which is also more reliable (75% of the unique hits come from LB; there is less variance). In addition, our results suggest, the effect of pH on synergy is weak in PA14.

We did Gene Set Enrichment Analysis (GSEA) using GO terms from pseudomonas.com website and Kologorov-Smirnov testing for statistical significance estimation. Although most common approach, it has been critizised for example here.

Cell component

library('topGO')
tab = getTopGO(chem.gen.res)

tmpFun = function (dat, name='CC') {
  dat[[name]][1:10, c('GO.ID', 'Term', 'raw.p.value')] %>% 
    mutate(Pvalue = round(as.numeric(raw.p.value), 4))  %>% 
    dplyr::select(-raw.p.value) %>% 
    mutate(`GO.ID` = sub('GO:', '', `GO.ID`)) 
}

tmpFun(tab) %>%
  knitr::kable()

GO.ID	Term	Pvalue
0008076	voltage-gated potassium channel complex	0.0096
0016021	integral component of membrane	0.0132
0055052	ATP-binding cassette (ABC) transporter complex, substrate-binding subunit-containing	0.0193
0005694	chromosome	0.0362
0033573	high-affinity iron permease complex	0.0530
0009289	pilus	0.1160
0005839	proteasome core complex	0.1221
0005960	glycine cleavage complex	0.1396
0030257	type III protein secretion system complex	0.1701
0005615	extracellular space	0.1717

  # knitr::kable(format='latex', booktabs=T) %>% 
  # kableExtra::save_kable('output/fig/table_cell_component.png', font='Verdana')

Biological process

# biological processes
tmpFun(tab, 'BP') %>% 
  knitr::kable()

GO.ID	Term	Pvalue
0009236	cobalamin biosynthetic process	0.0086
0009116	nucleoside metabolic process	0.0090
1900727	osmoregulated periplasmic glucan biosynthetic process	0.0166
0019700	organic phosphonate catabolic process	0.0173
0070475	rRNA base methylation	0.0179
0044249	cellular biosynthetic process	0.0260
0019354	siroheme biosynthetic process	0.0298
0006437	tyrosyl-tRNA aminoacylation	0.0346
0044262	cellular carbohydrate metabolic process	0.0357
0019748	secondary metabolic process	0.0370

  # knitr::kable(format='latex', booktabs=T) %>% 
  # kableExtra::save_kable('output/fig/table_biological_process.png', font='Verdana')

Protein-protein interaction network

Protein-protein interaction (PPI) analysis using STRING database. There was no data on PA14, so we will use PAO1 data to build and analyse the network onto which we then map PA14 orthologs.

library(igraph)
library(stringr)

hits = addAnnotation(chem.gen.res[t.test.q.value < 0.05, .(colony)])
hit_list = unique(hits[PAO1.ortholog != '', PAO1.ortholog])

# bring in the interactome
string_11_pa = 'input/dat/raw/pao1_ppi_stringDB_287.protein.links.v11.0.txt' %>% 
  fread() %>% 
  .[, `:=` (
    combined_score = combined_score/1000,
    #strip / trim species ID
    protein1 = gsub("287.DR97_", "PA", protein1, fixed=TRUE),
    protein2 = gsub("287.DR97_", "PA", protein2, fixed=TRUE)
  )]

# create an igraph object
pa_ppi = graph.data.frame(string_11_pa, directed = FALSE)
# summary(pa_ppi)
V(pa_ppi)$degree = degree(pa_ppi)

pa_ppi = simplify(pa_ppi, remove.multiple=T, remove.loops=T, edge.attr.comb='mean')

# extract network for all the hits in the provided list
hit_ppi = induced.subgraph(pa_ppi, which(V(pa_ppi)$name %in% hit_list))

# store degree on the network (main hit network)
V(hit_ppi)$degree = degree(hit_ppi)

# Replace the PAO1 genes with corresponding PA14 ones
lookup = unique(hits[, .(PAO1.ortholog, gene=gene.name.to.show)], by='PAO1.ortholog')
V(hit_ppi)$name = left_join(data.table(PAO1.ortholog=V(hit_ppi)$name), lookup)[, gene]

# Communities ------------------------------
# Community detection based on edge betweenness (Newman-Girvan)
# pdf('output/pdf/graph.pdf', height=10, width=20)

par(mfrow=c(1, 1), mar=c(1,3,1,1), cex=0.6)

set.seed(1)
ceb = cluster_edge_betweenness(hit_ppi) 
V(hit_ppi)$community = ceb$membership
cols = ins(cbPalette, 'white', c(5, 8:9, 11:14))

#svg('output/fig/SFig_ppi_cluster.svg', width=20, height=20)
#dendPlot(ceb, mode="hclust", palette = myCols, cex=1.5, cex.axis=1.5)
plot(hit_ppi, vertex.color=cols[V(hit_ppi)$community],
  vertex.label.cex=2, edge.width=4, vertex.label.family="Helvetica",
  vertex.label.dist=0,  vertex.size=16,
  vertex.frame.color='black', vertex.label.color='black'
)

#dev.off()

The major graph communities, using (Newman-Girvan’s edge betweenness):

1. The light blue nodes are central and seem to be enriched in regulatory genes.
2. The grey nodes seem to be a signal transduction from membrane to the regulatory genes (in light blue).
3. Orange and pink nodes are mostly related to metabolism; those in orange have more membrane related terms than the pink nodes.

Screen Validation

We validate the sensitivity of identified candidate mutants in low throughput and in liquid LB medium at pH 7.4. Instead of factorial (i.e. checkerboard), we use a fixed ratio design Tallarida et al 1997.

dat = fread('input/dat/raw/pa14_mut_pmb_rif.csv') %>% 
  addDrugRatio() %>% 
  addAUC() %>% 
  addFitness() %>% 
  remDissimilarExperimentalConditions() %>% 
  addPA14MutAnnotation()

dat[, dose:= ifelse(cond_f %in% c('PMB', 'combo'), d1, d2)]
dat_lst  = split(dat, dat$mut)
fit = lapply(dat_lst, fitPA14MutDR, curveid = cond_f)

# Loewe response surface ------------------------------
# rs = getPa14RS(dat)

# We'll use a snapshot of results from 2021-01-02
rsl = readRDS('input/dat/rds/2020-01-02_rs_pa14.rds') %>% lapply(., function(x) x$rsl)
foo = lapply(rsl, function(x) summary(x$maxR)$totals) %>% rbindlist(idcol='mut')
# Mutant numbering is potential source of confusion; hence I assign here
# manually: wild-type is coded as mut 0, fitting of mut 28 failed
foo[, mut := c(0:27, 29:45)]  
#arrange(foo, Syn) %>% print(20)

pltNSynergiesOfValidatedMutants()

# ggsave('output/fig/validated_mutants_synergy_count.svg', w=5, h=10)

Dose-response

#svg('output/fig/SFig_45PA14MutDoseResponses.svg', 12, 16)
#plt45PA14MutantDR()  # dose-responses
#dev.off()

Compare to Loewe’s null

#svg('output/fig/SFig_45PA14MutCompLoeweNull.svg', 12, 16)
# We'll use a snapshot of results from 2021-01-02
#plt45PA14MutantComparisonToLoeweNull('2020-01-02_rs_pa14')
#dev.off()

Table with E. coli orthologs

getValidHitTable() %>% addEcOrthologs() %>% .[,-'locus_pa'] %>% 
  knitr::kable()

PA gene	EC ortholog	Synergy	Description
02150		none	Serine phosphatase RsbU, regulator of sigma subunit
relA	relA	none	(p)ppGpp synthetase
26590		none	GntR family transcriptional regulator
66480		less	Predicted ATPase involved in chromosome partitioning
56840		less	Predicted thiol oxidoreductase
lldP	glcA	less	L-lactate permease
60490		less	cytochrome c
03760		more	sodium:solute symporter
43270	selU	more	tRNA 2-selenouridine synthase
51310		more	Predicted redox protein, regulator of disulfide bond formation
maiA		more	maleylacetoacetate isomerase
43670		more	sensor/response regulator hybrid
pmbA	pmbA	more	PmbA protein
62230		more	Predicted kinase

  # knitr::kable(format='latex', booktabs=T) %>% 
  # kableExtra::save_kable('output/fig/table_valid_hit.png')

GO term table

getGoTermTable() %>% 
  knitr::kable() 

PA gene	Location	Process	Function
02150	membrane	signal transduction	catalytic
03760	membrane	transmembrane transport	transmembrane transporter
26590	NA	regulation of transcription; biosynthesis	catalytic; transcription factor
43270	NA	tRNA seleno-modification	transferase for selenium-containing groups
43670	membrane	signal transduction; phosphorylation	phosphorelay sensor kinase
56840	NA	NA	electron transfer; heme binding
60490	NA	NA	electron transfer; heme binding
lldP	membrane	lactate transport	lactate transmembrane transport
maiA	cytoplasm	aromatic amino acid metabolism	catalytic; protein binding
pmbA	NA	peptidoglycan biosynthesis; proteolysis	metallopeptidase
relA	NA	ppGpp metabolic process	NA

  #knitr::kable(format='latex', booktabs=T) %>% 
  #kableExtra::save_kable('output/fig/table_valid_hit_go.png', font='Verdana')

Selected five mutants

3D

#rs = fread('input/dat/raw/pa14_mut_pmb_rif_cfu.csv') %>% 
#  prepForRS()  %>% 
#  setDT()  %>% 
#  addMutGeneNames() %>% 
#  getRS(., 'rs_6mut')
rs = readRDS('input/dat/tmp/rs_6mut.rds')

rsl = lapply(rs, function(x) lapply(x, function(y) y$rsl))

# lapply(names(rs), sav3D, xcap=10)

Membrane permeability

memb = readRDS('input/dat/rds/permeability.rds')
pltMemb(group_by(memb, PmbConcMgL, Strain))

#ggsave('./output/fig/PA14MembranePermeability.svg', w=15, h=5)

memb = readRDS('input/dat/rds/permeability.rds')
om = addAovDunnetMembrane(memb) 
im = addAovDunnetMembrane(memb, membrane = 'IM') 
getDunnetResultsMemb(om, lev = 1)

## $PmbConcMgL
## [1] 0.1
## 
## $RifConcMgL
## [1] 0
## 
## $Dunnett
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`wild-type`
##                      diff     lwr.ci    upr.ci   pval    
## 66480-wild-type -2.106667 -6.2881962 2.0748629 0.4809    
## 02150-wild-type -3.770000 -7.9515296 0.4115296 0.0819 .  
## relA-wild-type  -1.123333 -5.3048629 3.0581962 0.8942    
## 26590-wild-type -1.423333 -5.6048629 2.7581962 0.7821    
## 43270-wild-type  3.376667 -0.8048629 7.5581962 0.1304    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## $RobustDunnett
##  contrast            estimate   SE df t.ratio p.value
##  66480 - (wild-type)    -2.11 1.92 12  -1.098  0.6889
##  02150 - (wild-type)    -3.77 1.91 12  -1.974  0.2417
##  relA - (wild-type)     -1.12 3.01 12  -0.373  0.9773
##  26590 - (wild-type)    -1.42 1.92 12  -0.743  0.8698
##  43270 - (wild-type)     3.38 1.92 12   1.758  0.3293
## 
## P value adjustment: dunnettx method for 5 tests

getDunnetResultsMemb(im, lev = 1)

## $PmbConcMgL
## [1] 2
## 
## $RifConcMgL
## [1] 0
## 
## $Dunnett
## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $`wild-type`
##                       diff     lwr.ci      upr.ci    pval    
## 66480-wild-type  -7.536667 -16.605398   1.5320642  0.1162    
## 02150-wild-type  -9.753333 -18.822064  -0.6846025  0.0341 *  
## relA-wild-type  -25.990000 -35.058731 -16.9212692 2.2e-05 ***
## 26590-wild-type -29.013333 -38.082064 -19.9446025 5.3e-07 ***
## 43270-wild-type  10.973333   1.904602  20.0420642  0.0169 *  
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## $RobustDunnett
##  contrast            estimate   SE df t.ratio p.value
##  66480 - (wild-type)    -7.54 2.79 12  -2.704  0.0735
##  02150 - (wild-type)    -9.75 5.71 12  -1.708  0.3523
##  relA - (wild-type)    -25.99 3.07 12  -8.471  <.0001
##  26590 - (wild-type)   -29.01 2.14 12 -13.588  <.0001
##  43270 - (wild-type)    10.97 2.72 12   4.041  0.0069
## 
## P value adjustment: dunnettx method for 5 tests

Supplementary Membrane Assays

PMB resistant isolates

dat = 'input/dat/raw/pmb_resistant_isolates.csv' %>%  
  fread() %>% 
  setnames(., 'strain', 'mut') %>% 
  prepForRS(swapnames=T) %>% 
  setnames(., 'mut', 'strain') %>%
  setDT()

# rs = getRS(dat, 'rs_pmbR_intra')
# lapply(names(rs), sav3D)

Change in EC90 of RIF

dat = loadDdiData('input/dat/rds/ddi_all_strains.rds')

dat_fit = anti_join(dat, dat[c1==0 & c2!=0])[, .(data = list(.SD)), .(strain, cond)]

dat_fit[, fit := 
  map(data, function(x) {try(
    drm(effect ~ d1, d2, data = x, 
    fct=LL2.4(names=c("h", "m", "b", "e"), fixed=c(1, NA, NA, NA)),
    lowerl = c(-4.5, NA, NA))
  )})
]

## Error in optim(startVec, opfct, hessian = TRUE, method = "L-BFGS-B", lower = lowerLimits,  : 
##   non-finite value supplied by optim
## Error in drmOpt(opfct, opdfct1, startVecSc, optMethod, constrained, warnVal,  : 
##   Convergence failed

idx = sapply(dat_fit$fit, function(x) class(x) != 'try-error')

dat_fit[idx, ed90  := map(fit, function(x) list(try(ED(x, 90, interval='fls', display=F) %>% 
      round(3))))]

dat_fit[idx, par90 := lapply(unlist(dat_fit[['ed90']], recursive=F), function(x) try(as.data.table(x, keep.rownames='level')))]

par90 = lapply(1:nrow(dat_fit[idx, ]), 
  function(x) dat_fit[idx, ][['par90']][[x]]) %>% 
  rbindlist(idcol='id') %>% 
  merge(dat_fit[, .(strain, cond)][idx][, id := 1:26], .) %>% 
  lowerCaseNamesDT() %>% 
  .[, level := gsub('e:|:90', '', level)] %>% 
  .[, effect:= 90]

set.seed(1)
res = par90 %>% 
  group_by(id) %>% 
  mutate(rel_IC = estimate/head(estimate, 1),
    rel_IC = ifelse(rel_IC > 2, 2, rel_IC)
    ) %>% 
  mutate(level = ifelse(level < 1 & level != 0, 0.1, level)) %>% 
  filter(as.numeric(level) <= 1) 

# res %>% group_by(level) %>% summarise(ec90=median(estimate), quantile(estimate, 0.2), ci=gci(estimate), lwr=ec90/ci, upr=ec90*ci) %>% dplyr::select(-ci)
#   level  ec90   lwr   upr
# 1 0      7.41 6.34   8.65
# 2 0.1    3.20 2.32   4.41
# 3 1      2.18 0.966  4.92
# res %>% filter(level==0.1) %>% pull(estimate) %>% quantile(0.2) 
#   20% 
# 0.939 

ggplot(res, aes(level, estimate)) + 
  geom_jitter(width=0.2, pch=21, aes(fill=cond), col='grey70') +
  geom_boxplot(alpha=0.5, outlier.shape=NA) +
  scale_y_continuous(name=expression('EC'[90]^'RIF'*', mg/L'), 
    trans='log2', breaks=c(16, 8, 4, 2, 1, 0.5, 0.25)) +
  scale_x_discrete(name='PMB, xMIC', labels=c('0', '< 1', '1')) +
  coord_cartesian(ylim = c(0.2, 16)) +
  labs(fill = '') +
  theme(aspect.ratio=1)

#ggsave('output/fig/rif_ec90_upon_pmb.svg', width=5, height=4, bg='white')